﻿using System;
using System.Collections.Generic;
using System.Text;

namespace VerticalProcesses
{
    /// <summary>
    /// All radiation related equations go into here
    /// In case something needs to be calculated in the model
    /// the user can create this calculation by adding a series of equations
    /// from this folder and add his own ones into here
    /// (i.e. the user should break down long calculactions
    /// into single equations that go here)
    /// </summary>
    public static class RadiationEq
      
    {

        public static double ClearSkySolarRadiationOld(DateTime TimeStep, double TimeStepHrs, double latitude, double longitude, double elevation, double timezone, double slope, double aspect)
        {
            // Time dependent Variables
            //  hour = hour of the day
            //  min = minute
            //  dn = day number  (1-365)
            //  lapse = temp lapse rate (K/m)

            // Site characteristics
            //  elev = elevation (m)
            //  Latitude = Latitude (degree)
            //  Latitude 	 = 49 + (20/60) + (40/3600)
            //  Longitude = longitude (degree)    For convention W is negative, E is positive
            //  Longitude	 = -115 + (44/60) + (56/3600)
            //  Lstandard = standard longitude (degree)
            //  Lstandard = -120
            //  aspect = aspect at site   (degree)
            //  slope = slope inclination at site (degree)

            // Constants
            // I0 = solar constant (W/m^2)
            // I0 = 1368
            // pa = mean atmospheric clear sky transmissivity
            // pa = .75
            // g = 9.80665 (gravity - m/s^2)
            // M = 0.0289644(molar mass of earth's air kg/mol)
            // R = 8.31432 (universal gas constant N*m/mol*K)
            // Tb = 288.15 (Standard temp K)

            // CALIBRATION COEFFICIENTS
            // alphac = extinction coefficient for canopy
            // Stand Characteristics
            //dia = tree diameter (m)
            //hs = stem height (m)  
            //hc = height of the canopy (m)
            //n= tree density (stems/(m^2))
            double G0 = 0.0;
            double hour = 0.0;

            for (int i = 0; i < TimeStepHrs; i++)
            {
                //average for every interval > 1 hour
                if (TimeStepHrs == 1) { hour = TimeStep.Hour; }
                else { hour = (double)i; }
                double min = 30.0;      //middle of time step
                double dn = TimeStep.DayOfYear;

                double cr = Math.PI / 180.0;                //convert inputs to radians
                double a = (180.0 - aspect) * cr;  //convert to south=0, east  = positive as required in Iqbal equations
                double beta = slope * cr;
                double latRad = latitude * cr;
                double timeZoneRad = timezone * cr;        //Long = Longitude*con
                double hourDec = hour + min / 60.0;         //change to decimal hour

                double gamma = 2.0 * Math.PI * dn / 365.25;   //day angle (rad)       (Iqbal 1983)            
                double dcl = (0.006918 - 0.399912 * Math.Cos(gamma) + 0.070257 * Math.Sin(gamma)    //declination (Iqbal 1983)
                    - 0.006758 * Math.Cos(2.0 * gamma) + 0.000907 * Math.Sin(2.0 * gamma)
                    - 0.002697 * Math.Cos(3 * gamma) + 0.00148 * Math.Sin(3 * gamma));
                double Et = (0.000075 + 0.001868 * Math.Cos(gamma) - 0.032077 * Math.Sin(gamma) //Et = equation of time (hour) (Iqbal 1983)
                    - 0.014615 * Math.Cos(2.0 * gamma) - 0.04089 * Math.Sin(2.0 * gamma)) * 229.18 / 60;
                double latime = hourDec + (timezone - longitude) / 15.0 + Et;    //local apparent time (hour) (Iqbal 1983) 15 results from 4/60....

                //hour angle (noon zero, morning positive)      (from John Richards model 2008)
                double omega = (12.0 - latime) * Math.PI / 12.0;

                //solar zenith angle        (Iqbal 1983)
                double costheta = Math.Sin(dcl) * Math.Sin(latRad) + Math.Cos(dcl) * Math.Cos(latRad) * Math.Cos(omega);
                if (costheta < 0.0) { G0 = 0.0; continue; } //return 0.0 when sun position is below horizon
                double theta = Math.Acos(costheta);

                //solar azimuth angle (east positive, south = 0, west negative)(Iqbal 1983), costheta = sin solar altitude
                double cospsi = (costheta * Math.Sin(latRad) - Math.Sin(dcl)) / (Math.Sin(theta) * Math.Cos(latRad));
                double psi = Math.Acos(cospsi);
                if (omega < 0.0) { psi = 0.0 - psi; } //east positive, south = 0, west negative(Iqbal 1983)

                //cosine of angle of incidence
                double ci = Math.Cos(beta) * Math.Cos(theta) + Math.Sin(beta) * Math.Sin(theta) * Math.Cos(psi - a); //Eq 1.6.5b Iqbal 1983 == Eq 9 Moore and Rowland 1992
                if (ci < 0.0) { G0 = 0.0; continue; } //avoid 0 division

                //E0 = earth orbit eccentricity correction factor
                double E0 = 1.0 + 0.033 * Math.Cos(gamma - 0.048869);   //Iqbal 1983 and Grass Gis sun.r routine

                //barometric formula pp0=p/p0
                //double pp0 = Math.Pow((Constants.Tb / (Constants.Tb + tLapse * _XYZ[2])), ((Constants.g * Constants.M) / (Constants.R * tLapse)));
                //simpler formulation from r.sun in Grass:
                double pp0 = Math.Exp(-elevation / 8434.5);
                //k0 = potential direct open site radiation (W/m^2) on slope
                G0 += Constants.I0 * E0 * Math.Pow(Constants.pa, (pp0 / Math.Cos(theta))) * ci;  //Hock 1999
        
            }
            return Math.Max(G0 / TimeStepHrs, 0.0);
          }
        public static double CalcSolarAzimuthAngle(double latRad, double declination, double theta, double hourAngle)
        {
            double cospsi = (Math.Cos(theta) * Math.Sin(latRad) - Math.Sin(declination)) / (Math.Sin(theta) * Math.Cos(latRad));
            double psi = Math.Acos(cospsi);
            if (hourAngle < 0.0) { psi = 0.0 - psi; } //east positive, south = 0, west negative(Iqbal 1983)
            return psi;
        }
        /// <summary>
        /// Calculates Solar Zenith angle
        /// </summary>
        /// <param name="latRad"></param>
        /// <param name="dcl"></param>
        /// <param name="omega"></param>
        /// <returns>solar zenith angle</returns>
        public static double CalcSolarZenithAngle(double latRad, double dcl, double omega)
        {
            double costheta = Math.Sin(dcl) * Math.Sin(latRad) + Math.Cos(dcl) * Math.Cos(latRad) * Math.Cos(omega);
            return Math.Acos(costheta);
        }

        public static double ClearSkySolarRadiation(DateTime TimeStep, double TimeStepHrs, double latitude, double longitude, double elevation, double timezone, double slope, double aspect, double dbh, double nStems, double hStem, double hCrown, double ALPHAC)
          {
              // Time dependent Variables
              //  hour = hour of the day
              //  min = minute
              //  dn = day number  (1-365)
              //  lapse = temp lapse rate (K/m)

              // Site characteristics
              //  elev = elevation (m)
              //  Latitude = Latitude (degree)
              //  Latitude 	 = 49 + (20/60) + (40/3600)
              //  Longitude = longitude (degree)    For convention W is negative, E is positive
              //  Longitude	 = -115 + (44/60) + (56/3600)
              //  Lstandard = standard longitude (degree)
              //  Lstandard = -120
              //  aspect = aspect at site   (degree)
              //  slope = slope inclination at site (degree)

              // Constants
              // I0 = solar constant (W/m^2)
              // I0 = 1368
              // pa = mean atmospheric clear sky transmissivity
              // pa = .75
              // g = 9.80665 (gravity - m/s^2)
              // M = 0.0289644(molar mass of earth's air kg/mol)
              // R = 8.31432 (universal gas constant N*m/mol*K)
              // Tb = 288.15 (Standard temp K)

              // CALIBRATION COEFFICIENTS
              // alphac = extinction coefficient for canopy
              // Stand Characteristics
              //dia = tree diameter (m)
              //hs = stem height (m)  
              //hc = height of the canopy (m)
              //n= tree density (stems/(m^2))
              double Gf = 0.0;
              double hour = 0.0;

              for (int i = 0; i < TimeStepHrs; i++)
              {
                  //average for every interval > 1 hour
                  if (TimeStepHrs == 1) { hour = TimeStep.Hour; }
                  else { hour = (double)i; }
                  double min = 30.0;      //middle of time step
                  double dn = TimeStep.DayOfYear;

                  double cr = Math.PI / 180.0;                //convert inputs to radians
                  double a = (180.0 - aspect) * cr;  //convert to south=0, east  = positive as required in Iqbal equations
                  double beta = slope * cr;
                  double latRad = latitude * cr;
                  double timeZoneRad = timezone * cr;        //Long = Longitude*con
                  double hourDec = hour + min / 60.0;         //change to decimal hour

                  double gamma = 2.0 * Math.PI * dn / 365.25;   //day angle (rad)       (Iqbal 1983)            
                  double dcl = (0.006918 - 0.399912 * Math.Cos(gamma) + 0.070257 * Math.Sin(gamma)    //declination (Iqbal 1983)
                      - 0.006758 * Math.Cos(2.0 * gamma) + 0.000907 * Math.Sin(2.0 * gamma)
                      - 0.002697 * Math.Cos(3 * gamma) + 0.00148 * Math.Sin(3 * gamma));
                  double Et = (0.000075 + 0.001868 * Math.Cos(gamma) - 0.032077 * Math.Sin(gamma) //Et = equation of time (hour) (Iqbal 1983)
                      - 0.014615 * Math.Cos(2.0 * gamma) - 0.04089 * Math.Sin(2.0 * gamma)) * 229.18 / 60;
                  double latime = hourDec + (timezone - longitude) / 15.0 + Et;    //local apparent time (hour) (Iqbal 1983) 15 results from 4/60....

                  //hour angle (noon zero, morning positive)      (from John Richards model 2008)
                  double omega = (12.0 - latime) * Math.PI / 12.0;

                  //solar zenith angle        (Iqbal 1983)
                  double costheta = Math.Sin(dcl) * Math.Sin(latRad) + Math.Cos(dcl) * Math.Cos(latRad) * Math.Cos(omega);
                  if (costheta < 0.0) { Gf += 0.0; continue; } //return 0.0 when sun position is below horizon
                  double theta = Math.Acos(costheta);

                  //solar azimuth angle (east positive, south = 0, west negative)(Iqbal 1983), costheta = sin solar altitude
                  double cospsi = (costheta * Math.Sin(latRad) - Math.Sin(dcl)) / (Math.Sin(theta) * Math.Cos(latRad));
                  double psi = Math.Acos(cospsi);
                  if (omega < 0.0) { psi = 0.0 - psi; } //east positive, south = 0, west negative(Iqbal 1983)

                  //cosine of angle of incidence
                  double ci = Math.Cos(beta) * Math.Cos(theta) + Math.Sin(beta) * Math.Sin(theta) * Math.Cos(psi - a); //Eq 1.6.5b Iqbal 1983 == Eq 9 Moore and Rowland 1992
                  if (ci < 0.0) { Gf += 0.0; continue; } //avoid 0 division

                  //As = area of shadow from a stem  (m^2)
                  double As = dbh * hStem * Math.Sin(theta) / ci;     //Moore and Rowland 1992 Eq 8

                  //Lc = path length through canopy (m)
                  double Lc = hCrown * Math.Cos(beta) / ci; //Moore and Rowland 1992 Eq 7


                  //double tauCanopy = Math.Exp(-_alphac * _CrownClosure / Math.Cos(theta));

                  //E0 = earth orbit eccentricity correction factor
                  double E0 = 1.0 + 0.033 * Math.Cos(gamma - 0.048869);   //Iqbal 1983 and Grass Gis sun.r routine

                  //barometric formula pp0=p/p0
                  //double pp0 = Math.Pow((Constants.Tb / (Constants.Tb + tLapse * _XYZ[2])), ((Constants.g * Constants.M) / (Constants.R * tLapse)));
                  //simpler formulation from r.sun in Grass:
                  double pp0 = Math.Exp(-elevation / 8434.5);
                  //k0 = potential direct open site radiation (W/m^2) on slope
                  double G0 = Constants.I0 * E0 * Math.Pow(Constants.pa, (pp0 / Math.Cos(theta))) * ci;  //Hock 1999

                  //georg: Kf = potential k under forest (W/m^2), becomes kO at clearcuts since exp(0)=1
                  Gf += G0 * Math.Exp(-ALPHAC * Lc) * Math.Exp(-nStems * As); //Dan Moore personal communication
              
              }
              return Math.Max(Gf / TimeStepHrs, 0.0);
              
          }
        
        /// <summary>
         /// calculates the declination of the sun for a given julian day in rad
        /// </summary>
        /// <param name="julDay"></param>
        /// <returns>declination [rad]</returns>
          public static double CalcSunDeclination(int julDay)
          {
              double declination = 0.40954 * Math.Sin(0.0172 * (julDay - 79.35));
              return declination;
          }
        
        /// <summary>
          /// Calculates the solar constant sun for a given julian day
        /// </summary>
        /// <param name="julDay"></param>
        /// <returns>Solar constant [MJ/m2/min]</returns>
          public static double CalcSolarConstant(int julDay)
          {
              //solar constant in J / m² / min
              double S = 81930 + 2910 * Math.Cos(Math.PI / 180 * (julDay - 15));
              // J --> MJ
              S = S / 1000000;

              return S;
          }
        
        /// <summary>
        /// calculates the inverse relative distance Earth-Sun in rad
        /// </summary>
        /// <param name="julDay"></param>
          /// <returns>inverse relative distance Earth-sun [rad]</returns>
          public static double CalcInverseRelativeDistanceEarthSun(int julDay)
          {
              double irdes = 1 + 0.033 * Math.Cos((2 * Math.PI / 365) * julDay);
              return irdes;
          }

        /// <summary>
        /// calculates the sunset hour angle in rad
        /// </summary>
        /// <param name="latitude"></param>
        /// <param name="declination"></param>
          /// <returns></returns>
          public static double CalcSunsetHourAngleDaily(double latitude, double declination)
          {
              double sha = Math.Acos(-1 * Math.Tan(latitude) * Math.Tan(declination));
              return sha;
          }

        /// <summary>
          /// calculates the mid time hour angle between two time steps
        /// </summary>
        /// <param name="hour"></param>
        /// <param name="julDay"></param>
          /// <param name="longSite">longitude of the point of interest in deg</param>
        /// <param name="longTZ"></param>
          /// <returns>hour angle in [deg]</returns>
        public static double CalcMidTimeHourAngle(long hour, int julDay, double longSite, double longTZ){
        long minute = 0;
        double decMin = minute / 60;
        double decTime = hour + decMin;
        if(decTime >= 24.0)
            decTime = decTime - 24;
        double midTime = decTime + 0.5;
        double b = (2 * Math.PI * (julDay - 81)) / 364;
        double Sc = 0.1645 * Math.Sin(2*b)-0.1255 * Math.Cos(b) - 0.025 * Math.Sin(b);
        double midTimeHourAngle = Math.PI / 12 * (12 - (midTime + 0.06667 * (longTZ - longSite) + Sc));//Georg: there was a mistake (i think) with -12...
        return midTimeHourAngle;
    }

        /// <summary>
        /// calculates the hour angle at the beginning of the current hourly time step
        /// </summary>
        /// <param name="midTimeHourAngle"></param>
        /// <returns></returns>
          public static double CalcStartTimeHourAngle(double midTimeHourAngle)
          {
              double startTimeHourAngle = (midTimeHourAngle - (Math.PI * 1 / 24));
              return startTimeHourAngle;
          }

        /// <summary>
        /// calculates the hour angle at the end of the current hourly time step
        /// </summary>
        /// <param name="midTimeHourAngle"></param>
        /// <returns></returns>
          public static double CalcEndTimeHourAngle(double midTimeHourAngle)
          {
              double endTimeHourAngle = (midTimeHourAngle + (Math.PI * 1 / 24));
              return endTimeHourAngle;
          }


        /// <summary>
        /// calculates the daily extraterrestrial radiation
        /// </summary>
        /// <param name="Gsc">solar constant [MJ / m²/ min]</param>
        /// <param name="dr">inverse relative distance Earth-Sun [rad.]</param>
        /// <param name="ws">hour angle [rad.]</param>
        /// <param name="lat">latitude of the point of interest [rad.]</param>
        /// <param name="decl">sun's declination [rad.</param>
        /// <returns>extraterrestrial radiation [MJ / m² day]</returns>
          public static double CalcDailyExtraterrestrialRadiation(double Gsc, double dr,double ws, double lat, double decl)
          {
              double Ra = ((24 * 60) / Math.PI) * Gsc * dr * (ws * Math.Sin(lat) * Math.Sin(decl) + Math.Cos(lat) * Math.Cos(decl) * Math.Sin(ws));
              return Ra;
          }
        /// <summary>
        /// calculates the hourly extraterrestrial radiation
        /// </summary>
        /// <param name="Gsc">solar constant of the current time step [MJ/m2/min]</param>
        /// <param name="relDist">relative distance Earth -- Sun [rad]</param>
        /// <param name="w1">hour angle at the beginning of the time step [rad]</param>
        /// <param name="w2">hour angle at the end of the time step [rad]</param>
        /// <param name="radLat"></param>
        /// <param name="decl"></param>
        /// <returns>extraterrestrial radiation [MJ / m²/ hour]</returns>
          public static double CalcHourlyExtraterrestrialRadiation(double Gsc, double relDist, double w1, double w2, double radLat, double decl)
          {
              double Ra = ((12 * 60) / Math.PI) * Gsc * relDist * ((w2 - w1) * Math.Sin(radLat) * Math.Sin(decl) + Math.Cos(radLat) * Math.Cos(decl) * (Math.Sin(w2) - Math.Sin(w1)));
              //no negative radiation !!
              if (Ra < 0)
                  Ra = 0;
              return Ra;
          }
     /// <summary>
          /// calculates the clear sky solar radiation in MJ /m²/timestep (daily or hourly)
     /// <param name="timeStep"></param>
     /// <param name="TimeStepHrs"></param>
     /// <param name="latitude"></param>
     /// <param name="longitude"></param>
     /// <param name="elevation"></param>
     /// <param name="timezone"></param>
     /// <param name="slope"></param>
     /// <param name="aspect"></param>
          /// <returns>clear sky solar radiation MJ /m²/timestep</returns>
          public static double CalcClearSkySolarRadiation(DateTime timeStep, double timeStepHrs, double latitude, double longitude, double elevation, double timezone, double slope, double aspect)
          {
              double extraTerrRad = 0.0;
              int hour = timeStep.Hour;
              int julDay = timeStep.DayOfYear;
              //double conv = 1.0 / (3600.0 * timeStepHrs) * 1.0E6; //convert MJ to W  (60 [sec/min] * 60 [min/h])
              double latRad = latitude * (Math.PI / 180.0);
              double declination = CalcSunDeclination(julDay);
              double solarConstant = CalcSolarConstant(julDay);
              double invRelDistEarthSun = CalcInverseRelativeDistanceEarthSun(julDay);
              double SACF = CalcSlopeAspectCorrectionFactor(latRad, slope, aspect, declination);
              if (timeStepHrs == 24)
              {
                  double sunsetHourAngle = CalcSunsetHourAngleDaily(latRad, declination);
                  extraTerrRad = CalcDailyExtraterrestrialRadiation(solarConstant, invRelDistEarthSun, sunsetHourAngle, latRad, declination);
                  //for (int i=0; i < 24; i++ )
                  //{
                  //    double midTimeHourAngle = CalcMidTimeHourAngle(i, julDay, longitude, timezone);
                  //    double startTimeHourAngle = CalcStartTimeHourAngle(midTimeHourAngle);
                  //    double endTimeHourAngle = CalcEndTimeHourAngle(midTimeHourAngle);
                  //    extraTerrRad += CalcHourlyExtraterrestrialRadiation(solarConstant, invRelDistEarthSun, startTimeHourAngle, endTimeHourAngle, latRad, declination);
                  //    }
              }
              else //hourly
              {
                  double midTimeHourAngle = CalcMidTimeHourAngle(hour, julDay, longitude, timezone);
                  double startTimeHourAngle = CalcStartTimeHourAngle(midTimeHourAngle);
                  double endTimeHourAngle = CalcEndTimeHourAngle(midTimeHourAngle);
                  extraTerrRad = CalcHourlyExtraterrestrialRadiation(solarConstant, invRelDistEarthSun, startTimeHourAngle, endTimeHourAngle, latRad, declination);
              }

              return extraTerrRad * SACF * (0.75 + 2E-5 * elevation);

          }
        /// <summary>
          /// Extraterrestrial Radiation in MJ/timestep (hourly or daily)
        /// </summary>
        /// <param name="timeStep">current date and time</param>
        /// <param name="TimeStepHrs">hours/timestep</param>
          /// <param name="latitude">latitude</param>
        /// <param name="longitude"></param>
        /// <param name="elevation"></param>
        /// <param name="timezone"></param>
        /// <param name="slope"></param>
        /// <param name="aspect"></param>
          /// <returns>Extraterrestrial Radiation in MJ/timestep (hourly or daily)</returns>
        public static double CalcExtraterrRad(DateTime timeStep,double TimeStepHrs, double latitude, double longitude, double elevation, double timezone, double slope, double aspect)
        {
            double extraTerrRad = 0.0;
            int hour = timeStep.Hour;
            int julDay = timeStep.DayOfYear;
            double latRad = latitude * (Math.PI / 180.0);
            double declination = CalcSunDeclination(julDay);
            double solarConstant = CalcSolarConstant(julDay);
            double invRelDistEarthSun = CalcInverseRelativeDistanceEarthSun(julDay);
            double SACF = CalcSlopeAspectCorrectionFactor(latRad, slope, aspect, declination);
            if (TimeStepHrs == 24)
            {
                double sunsetHourAngle = CalcSunsetHourAngleDaily(latRad, declination);
                extraTerrRad = CalcDailyExtraterrestrialRadiation(solarConstant, invRelDistEarthSun, sunsetHourAngle, latRad, declination);
            }
            else //hourly
            {
                double midTimeHourAngle = CalcMidTimeHourAngle(hour, julDay, longitude, timezone);
                double startTimeHourAngle = CalcStartTimeHourAngle(midTimeHourAngle);
                double endTimeHourAngle = CalcEndTimeHourAngle(midTimeHourAngle);
                int idx = julDay * 24 + hour;
                extraTerrRad = CalcHourlyExtraterrestrialRadiation(solarConstant, invRelDistEarthSun, startTimeHourAngle, endTimeHourAngle, latRad, declination);
            }
            return extraTerrRad * SACF;

        }

       
       /// <summary>
        /// calculates the influence of the slope-aspect combination of point of interest
       /// </summary>
        /// 
       /// <param name="latRad"></param>
       /// <param name="slope"></param>
       /// <param name="aspect"></param>
       /// <param name="declRad"></param>
       /// <returns>slope aspect corretion factor</returns>
        public static double CalcSlopeAspectCorrectionFactor(double latRad, double slope, double aspect, double declRad)
        {
            double convAsp = 180 - aspect;
            double aspRad = convAsp * Math.PI / 180.0;
            double slopeRad = slope * Math.PI / 180.0;

            double sloped = (Math.Sin(declRad) * Math.Sin(latRad) * Math.Cos(slopeRad))
            - (Math.Sin(declRad) * Math.Cos(latRad) * Math.Sin(slopeRad) * Math.Cos(aspRad))
            + (Math.Cos(declRad) * Math.Cos(latRad) * Math.Cos(slopeRad) * 1)
            + (Math.Cos(declRad) * Math.Sin(latRad) * Math.Sin(slopeRad) * Math.Cos(aspRad) * 1)
            + (Math.Cos(declRad) * Math.Sin(slopeRad) * Math.Sin(aspRad) * 0);
            double horizontal = (Math.Sin(declRad) * Math.Sin(latRad) + Math.Cos(declRad) * Math.Cos(latRad) * 1);
            double slopeAspect = sloped / horizontal;

            if (slopeAspect < 0) //can happen
                slopeAspect = 0;

            return slopeAspect;
        }

/// <summary>
/// Calculates radation below canopy according to Rowland and Moore
/// solar angles are calcualted according to Iqbal
/// notation of Iqbal and Rowland and Moore
/// </summary>
/// <param name="TimeStep"></param>
/// <param name="Radiation"></param>
/// <param name="TimeStepHrs"></param>
/// <param name="latitude"></param>
/// <param name="longitude"></param>
/// <param name="elevation">m</param>
/// <param name="timezone"></param>
/// <param name="slope">degrees</param>
/// <param name="aspect">degrees</param>
/// <param name="dbh">m</param>
/// <param name="nStems">stems/m2</param>
/// <param name="hStem">m</param>
/// <param name="hCrown">m</param>
/// <param name="ALPHAC"></param>
/// <returns>Radiation below canopy [W/m2]</returns>
        public static double CalcRadBelowCanopy(DateTime TimeStep, double Radiation, double TimeStepHrs, double latitude, double longitude, double elevation, 
            double timezone, double slope, double aspect, double dbh, double nStems, double hStem, double hCrown, double ALPHAC)
        {
            double Kf = 0.0;
            int hour = TimeStep.Hour;
            int julDay = TimeStep.DayOfYear;
            double conv = (Math.PI / 180.0);
            double latRad = latitude * conv;
            double slopeRad = slope * conv;
            double aspectRad = aspect * conv;
            double declination = CalcSunDeclination(julDay);           
            if (TimeStepHrs == 24)
            {
                double midTimeHourAngle = CalcMidTimeHourAngle(12, julDay, longitude, timezone);
                double theta = CalcSolarZenithAngle(latRad, declination, midTimeHourAngle);
                double psi = CalcSolarAzimuthAngle(latRad, declination, theta, midTimeHourAngle);
                double ci = CalcSolarIncidenceAngle(aspectRad, slopeRad, theta, psi);
                double As = dbh * hStem * Math.Sin(theta) / ci;
                double lc = hCrown * Math.Cos(slopeRad) / ci;
                Kf = Radiation * Math.Exp(-ALPHAC * lc) * Math.Exp(-nStems * As);
            }
            else //hourly
            {
                double midTimeHourAngle = CalcMidTimeHourAngle(hour, julDay, longitude, timezone);
                double theta = CalcSolarZenithAngle(latRad, declination, midTimeHourAngle);
                double psi = CalcSolarAzimuthAngle(latRad, declination, theta, midTimeHourAngle);
                double ci = CalcSolarIncidenceAngle(aspectRad, slopeRad, theta, psi);
                if (ci < 0.0) { return Kf = 0.0; }
                double As = dbh * hStem * Math.Sin(theta) / ci;
                double lc = hCrown * Math.Cos(slopeRad) / ci;
                Kf = Radiation * Math.Exp(-ALPHAC * lc) * Math.Exp(-nStems * As);
              }
            return Kf;

        }

        public static double CalcSolarIncidenceAngle(double aspectRad, double slopeRad, double theta, double psi)
        {
            double a = Math.PI - aspectRad; //convert to south=0, east  = positive as required in Iqbal equations
            return Math.Cos(slopeRad) * Math.Cos(theta) + Math.Sin(slopeRad) * Math.Sin(theta) * Math.Cos(psi - a);

        }
        /// <summary>
        /// Adjust measured radiation to slope and aspect for a given location
        /// </summary>
        /// <param name="timeStep"></param>
        /// <param name="radiation"></param>
        /// <param name="latitude"></param>
        /// <param name="slope"></param>
        /// <param name="aspect"></param>
        /// <returns>slopeAspect adjusted radiation [W/m2]</returns>
        public static double CalcRadSlopeAspect(DateTime timeStep, double radiation, double latitude, double slope, double aspect)
        {
            int julDay = timeStep.DayOfYear;
            double conv = (Math.PI / 180.0);
            double latRad = latitude * conv;
            double declination = CalcSunDeclination(julDay);
            double SACF = CalcSlopeAspectCorrectionFactor(latRad, slope, aspect, declination);
            return radiation * SACF;
        }
        /// <summary>
        /// calculates the daily net shortwave radiation resulting from the balance between incoming
        /// and reflected solar radiation
        /// </summary>
        /// <param name="albedo">albedo of the landcover</param>
        /// <param name="Rs">Rs the daily solar radiation [MJ / m² day or hour]</param>
        /// <returns>net solar or shortwave radiation [MJ / m² day or hour]</returns>
        public static double CalcNetShortwaveRadiation(double albedo, double Rs)
        {
            double Rns = (1 - albedo) * Rs;
            return Rns;
        }
        /// <summary>
        /// calculates the net radiation
        /// </summary>
        /// <param name="Rns">daily net solar or shortwave radiation [MJ / m² day or hour]</param>
        /// <param name="Rnl">Rnl the daily net longwave radiation [MJ / m² day or hour]</param>
        /// <returns>daily net radiation [MJ / m² day] or hour</returns>
        public static double CalcNetRadiation(double Rns, double Rnl)
        {
            double Rn = Rns - Rnl;
            if (Rn < 0)
                Rn = 0;
            return Rn;
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="tmean"></param>
        /// <param name="ea"></param>
        /// <param name="Rs"></param>
        /// <param name="Rs0"></param>
        /// <param name="debug"></param>
        /// <returns></returns>
        public static double CalcDailyNetLongwaveRadiation(double Tmax, double Tmin, double Ea, double Rs, double Rso)
        {
        /** the Stefan Bolztmann constant in [MJ / K^4 m² day] **/
            double Rnl = Constants.StefanBolzConstMJ * (Math.Pow((Tmax + 237.15), 4) + Math.Pow((Tmin + 237.16), 4)) * (0.34 - 0.14 * Math.Sqrt(Ea)) * (1.35 * Rs / Rso - 0.35) / 2;
        return Rnl;
    }
    }
}
